After viewing the size factors, it appears that we have 5 samples with really low values that should be removed:
Without filtering low reads, we had a total of 38828 counts. After filtering out the low counts (those with a base mean less than 3), we now have 24384 counts remaining for all C. virginica samples.
# going to use the filter_counts_out object created above because it is filtered for low reads and has outliers removed
## create the DESeq object
dds <- DESeqDataSetFromMatrix(countData = filter_counts_out,
colData = expDesign,
design = ~ treat)
dds <- DESeq(dds) # differential expression analysis on gamma-poisson distribution
vsd <-varianceStabilizingTransformation(dds, blind = TRUE) # quickly estimate dispersion trend and apply a variance stabilizing transformation
## Save dds and vsd data
save(dds, vsd, file = "Data/transformed_counts.RData")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
##
## adonis2(formula = pcaData[6:31] ~ pcaData$infect * pcaData$pCO2, permutations = 1500, method = "eu")
## Df SumOfSqs R2 F Pr(>F)
## pcaData$infect 1 1118.3 0.04736 1.1893 0.1799
## pcaData$pCO2 1 945.3 0.04003 1.0053 0.3817
## pcaData$infect:pcaData$pCO2 1 864.0 0.03659 0.9188 0.5623
## Residual 22 20686.8 0.87603
## Total 25 23614.3 1.00000
PCA plot with each point being one oyster sample colored by treatment. Orange represents the control treatment (no boring sponge, pCO2 level 400), green represents no sponge, pCO2 level 2800, purple represents samples infected with boring sponge at pCO2 2800, and pink represents samples infected with boring sponge at pCO2 400.
## log2 fold change (MLE): treat S 400 vs N 2800
## Wald test p-value: treat S 400 vs N 2800
## DataFrame with 24384 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## LOC111126949 3.452259 -0.9422247 0.817784 -1.1521686 0.249252 0.999652
## LOC111110729 0.485466 -0.5124481 1.463646 -0.3501175 0.726251 0.999652
## LOC111120752 1.693621 -0.2827063 0.846812 -0.3338478 0.738494 0.999652
## LOC111113860 5.933380 -0.0448271 0.601266 -0.0745544 0.940569 0.999652
## LOC111109550 0.331392 -0.9860539 2.483653 -0.3970176 0.691355 0.999652
## ... ... ... ... ... ... ...
## LOC111117112 0.6463025 1.6788746 1.80337 0.9309629 0.351873 0.999652
## LOC111117113 0.1436384 0.7681374 3.83879 0.2000987 0.841403 0.999652
## LOC111117117 0.0778208 0.2872407 4.24705 0.0676329 0.946078 0.999652
## LOC111116908 0.1097457 -0.0817811 4.24580 -0.0192617 0.984632 0.999652
## LOC111117715 0.6210268 -0.4455042 1.67997 -0.2651862 0.790866 0.999652
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 2, 0.0082%
## LFC < 0 (down) : 0, 0%
## outliers [1] : 36, 0.15%
## low counts [2] : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_con_sponge)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 0, 0%
## LFC < 0 (down) : 0, 0%
## outliers [1] : 36, 0.15%
## low counts [2] : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_con_pco2)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 0, 0%
## LFC < 0 (down) : 0, 0%
## outliers [1] : 36, 0.15%
## low counts [2] : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_trt)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 6, 0.025%
## LFC < 0 (down) : 5, 0.021%
## outliers [1] : 36, 0.15%
## low counts [2] : 11339, 47%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_trt_pco2)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 3, 0.012%
## LFC < 0 (down) : 5, 0.021%
## outliers [1] : 36, 0.15%
## low counts [2] : 14166, 58%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_stn_pco2)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 2, 0.0082%
## LFC < 0 (down) : 0, 0%
## outliers [1] : 36, 0.15%
## low counts [2] : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_shn_pco2)
##
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 1, 0.0041%
## LFC < 0 (down) : 3, 0.012%
## outliers [1] : 36, 0.15%
## low counts [2] : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
It is just a demo, not very exciting :D
Session information from the last full knit of Rmarkdown on 30 June 2022.
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggvenn_0.1.9 adegenet_2.1.7
## [3] ade4_1.7-19 ggrepel_0.9.1
## [5] pdftools_3.2.1 ggpubr_0.4.0
## [7] data.table_1.14.2 vegan_2.6-2
## [9] lattice_0.20-45 permute_0.9-7
## [11] DESeq2_1.36.0 SummarizedExperiment_1.26.1
## [13] Biobase_2.56.0 MatrixGenerics_1.8.1
## [15] matrixStats_0.62.0 GenomicRanges_1.48.0
## [17] GenomeInfoDb_1.32.2 IRanges_2.30.0
## [19] S4Vectors_0.34.0 BiocGenerics_0.42.0
## [21] plotly_4.10.0 arrayQualityMetrics_3.52.0
## [23] forcats_0.5.1 stringr_1.4.0
## [25] dplyr_1.0.9 purrr_0.3.4
## [27] readr_2.1.2 tidyr_1.2.0
## [29] tibble_3.1.7 ggplot2_3.3.6
## [31] tidyverse_1.3.1 knitr_1.39
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 RSQLite_2.2.14
## [4] AnnotationDbi_1.58.0 htmlwidgets_1.5.4 beadarray_2.46.0
## [7] BiocParallel_1.30.3 munsell_0.5.0 codetools_0.2-18
## [10] preprocessCore_1.58.0 withr_2.5.0 colorspace_2.0-3
## [13] highr_0.9 rstudioapi_0.13 setRNG_2022.4-1
## [16] ggsignif_0.6.3 labeling_0.4.2 GenomeInfoDbData_1.2.8
## [19] hwriter_1.3.2.1 bit64_4.0.5 farver_2.1.0
## [22] vctrs_0.4.1 generics_0.1.2 xfun_0.31
## [25] qpdf_1.2.0 R6_2.5.1 illuminaio_0.38.0
## [28] gridSVG_1.7-4 locfit_1.5-9.5 bitops_1.0-7
## [31] cachem_1.0.6 DelayedArray_0.22.0 assertthat_0.2.1
## [34] promises_1.2.0.1 scales_1.2.0 nnet_7.3-17
## [37] gtable_0.3.0 affy_1.74.0 rlang_1.0.2
## [40] genefilter_1.78.0 systemfonts_1.0.4 splines_4.2.0
## [43] rstatix_0.7.0 lazyeval_0.2.2 hexbin_1.28.2
## [46] broom_0.8.0 checkmate_2.1.0 BiocManager_1.30.18
## [49] yaml_2.3.5 reshape2_1.4.4 abind_1.4-5
## [52] modelr_0.1.8 crosstalk_1.2.0 backports_1.4.1
## [55] httpuv_1.6.5 Hmisc_4.7-0 tools_4.2.0
## [58] affyio_1.66.0 ellipsis_0.3.2 jquerylib_0.1.4
## [61] RColorBrewer_1.1-3 Rcpp_1.0.8.3 plyr_1.8.7
## [64] base64enc_0.1-3 zlibbioc_1.42.0 RCurl_1.98-1.7
## [67] rpart_4.1.16 openssl_2.0.2 haven_2.5.0
## [70] cluster_2.1.3 fs_1.5.2 magrittr_2.0.3
## [73] reprex_2.0.1 hms_1.1.1 mime_0.12
## [76] evaluate_0.15 xtable_1.8-4 XML_3.99-0.10
## [79] jpeg_0.1-9 gcrma_2.68.0 readxl_1.4.0
## [82] gridExtra_2.3 compiler_4.2.0 crayon_1.5.1
## [85] htmltools_0.5.2 mgcv_1.8-40 later_1.3.0
## [88] tzdb_0.3.0 Formula_1.2-4 geneplotter_1.74.0
## [91] lubridate_1.8.0 DBI_1.1.3 dbplyr_2.2.1
## [94] MASS_7.3-57 Matrix_1.4-1 car_3.1-0
## [97] cli_3.3.0 vsn_3.64.0 parallel_4.2.0
## [100] igraph_1.3.2 pkgconfig_2.0.3 foreign_0.8-82
## [103] xml2_1.3.3 svglite_2.1.0 annotate_1.74.0
## [106] bslib_0.3.1 BeadDataPackR_1.48.0 affyPLM_1.72.0
## [109] XVector_0.36.0 rvest_1.0.2 digest_0.6.29
## [112] Biostrings_2.64.0 rmarkdown_2.14 base64_2.0
## [115] cellranger_1.1.0 htmlTable_2.4.0 shiny_1.7.1
## [118] lifecycle_1.0.1 nlme_3.1-158 jsonlite_1.8.0
## [121] carData_3.0-5 seqinr_4.2-16 viridisLite_0.4.0
## [124] askpass_1.1 limma_3.52.2 fansi_1.0.3
## [127] pillar_1.7.0 KEGGREST_1.36.2 fastmap_1.1.0
## [130] httr_1.4.3 survival_3.3-1 glue_1.6.2
## [133] png_0.1-7 bit_4.0.4 stringi_1.7.6
## [136] sass_0.4.1 blob_1.2.3 latticeExtra_0.6-29
## [139] memoise_2.0.1 ape_5.6-2